function [model]=r_nemo_mx(dat,parameter,bbio,bbio_dir,year)  


switch parameter
    case 'T'
        gridtype='gridT';
        kind='T';
        levmin=1;
        if (length(fieldnames(dat)) < 5)  % no depth - surface
            levmax=1;
            [model.lon,model.lat,model.data,time_sec,time_txt,deptht]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax);
        else
            levmax=28; % Clim only to 1500m ...  then lev 28 enough
            [model.lon,model.lat,model.data,time_sec,time_txt,model.depth]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax,year);    
            model.depth=model.depth(1:levmax);
        end
        % time_txt -> matlab time
        model.time=datenum(str2num(char(time_txt(:,2:5))),str2num(char(time_txt(:,7:8))),str2num(char(time_txt(:,10:11))));
    case 'S'
        gridtype='gridT';
        kind='S';
        levmin=1;
        levmax=28; % Clim only to 1500m ...  then lev 28 enough
        [model.lon,model.lat,model.data,time_sec,time_txt,model.depth]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax,year);   
        model.depth=model.depth(1:levmax);
        
        % time_txt -> matlab time
        model.time=datenum(str2num(char(time_txt(:,2:5))),str2num(char(time_txt(:,7:8))),str2num(char(time_txt(:,10:11))));
        
    case 'NO3'
        gridtype='ptrcT';
        kind='no3';
        levmin=1;
        levmax=28; % Clim only to 1500m ...  then lev 28 enough
        [model.lon,model.lat,model.data,time_sec,time_txt,model.depth]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax,year);   
        model.depth=model.depth(1:levmax);
        
        % time_txt -> matlab time
        model.time=datenum(str2num(char(time_txt(:,2:5))),str2num(char(time_txt(:,7:8))),str2num(char(time_txt(:,10:11))));

    case 'CHL'
        gridtype='ptrcT';
        levmin=1;
        levmax=1;
        kind='phy';
        [model.lon,model.lat,phy,time_sec,time_txt,model.depth]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax,year); 
        kind='zle';
        [model.lon,model.lat,zle,time_sec,time_txt,model.depth]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax,year);  
        kind='zlno3'; 
        [model.lon,model.lat,zlno3,time_sec,time_txt,model.depth]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax,year);     
        model.depth=model.depth(1:levmax);
        
    
        % Phytoplankton (mmolN.m-3) to Chlorophyll (mgChl.m-3) - variable
        % ratio
        interm(1:size(zle,1),1:size(zle,2),1:size(zle,3),1:size(zle,4))=1.;
        index=find(zle > zlno3);
        interm(index)=zlno3(index)./zle(index);
        interm(find(interm<0.31))=0.31;
        model.data=squeeze(1.59.*interm.*phy);      
        
        % Constant ratio
        %model.data=1.59.*squeeze(phy);
        
        % time_txt -> matlab time
        model.time=datenum(str2num(char(time_txt(:,2:5))),str2num(char(time_txt(:,7:8))),str2num(char(time_txt(:,10:11))));
        
    case 'OWP'
        gridtype='ptrcT';
      
        % Commented: option to reduce needed memory
        chl=zeros([17 529 486 73]);
        for ilev=1:17
            levmin=ilev;
            levmax=ilev; % 0-519m
        %    levmin=1;
        %    levmax=17; % 0-519m
            kind='phy';
            [model.lon,model.lat,phy,time_sec,time_txt,model.depth]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax,year); 
            kind='zle';
            [model.lon,model.lat,zle,time_sec,time_txt,model.depth]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax,year);  
            kind='zlno3'; 
            [model.lon,model.lat,zlno3,time_sec,time_txt,model.depth]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax,year);     
            model.depth=model.depth(1:levmax);
            
            % Phytoplankton (mmolN.m-3) to Chlorophyll (mgChl.m-3) - variable
            % ratio
            interm=ones([size(zle,1) size(zle,2) size(zle,3) size(zle,4)]);
            index=find(zle > zlno3);
            interm(index)=zlno3(index)./zle(index);
            interm(find(interm<0.31))=0.31;
        %    chl=squeeze(1.59.*interm.*phy);
            chl(ilev,:,:,:)=squeeze(1.59.*interm.*phy);
        end

        % Compute the optically weighted pigment
        kw=1/25;

        int1(1:size(chl,2),1:size(chl,3),1:size(chl,4))=0.;
        int2(1:size(chl,2),1:size(chl,3),1:size(chl,4))=0.;

        int1(:,:,:)=squeeze(chl(1,:,:,:)).*exp(-2*kw*model.depth(1)).*model.depth(1);
        int2(:,:,:)=exp(-2*kw*model.depth(1))*model.depth(1);        


        for ilev=2:17
          int1(:,:,:)=int1(:,:,:)+(squeeze(chl(ilev,:,:,:)).*exp(-2*kw*model.depth(ilev)).*(model.depth(ilev)-model.depth(ilev-1)));
          int2(:,:,:)=int2(:,:,:)+(exp(-2*kw*model.depth(ilev))*(model.depth(ilev)-model.depth(ilev-1)));    
        end

        model.data(1:size(chl,2),1:size(chl,3),1:size(chl,4))=NaN;
        index=find(int2 > 0.);
        model.data(index)=int1(index)./int2(index);

        % time_txt -> matlab time
        model.time=datenum(str2num(char(time_txt(:,2:5))),str2num(char(time_txt(:,7:8))),str2num(char(time_txt(:,10:11))));
       
        
    case 'PP'
        gridtype='ptrcT';
        kind='PP';
        levmin=1;
        levmax=12;
        [model.lon,model.lat,pp,time_sec,time_txt,model.depth]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax,year);
        
        model.data=squeeze(pp(1,:,:,:)).*model.depth(1);
        for jpk=2:12
            model.data=model.data+(squeeze(pp(jpk,:,:,:)).*(model.depth(jpk)-model.depth(jpk-1)));       
        end
        %mmolN.m-2.s-1 -> mgC.m-2.j-1
        model.data=model.data.*2400.*6.625.*12./5;
        
        % time_txt -> matlab time
        model.time=datenum(str2num(char(time_txt(:,2:5))),str2num(char(time_txt(:,7:8))),str2num(char(time_txt(:,10:11))));

        
    case 'EKE'
        % ---- Read Model output
        
        % -- From U,V
        [lon,lat,u,time_sec,time_txt,deptht]=rnemo_3d(bbio,bbio_dir,'gridU','u',1,1);
        [lon,lat,v,time_sec,time_txt,deptht]=rnemo_3d(bbio,bbio_dir,'gridV','v',1,1);
        
        u=squeeze(u).*100.;
        v=squeeze(v).*100.;
        
        % ---- Compute the temporal mean and temporal quadratic mean:
        um=mean(u,3);
        vm=mean(v,3);
        u2m=mean(u.*u,3);
        v2m=mean(v.*v,3);
        
        % ---- Compute the EKE
        eke_tot(1:size(u,1),1:size(v,2))=0.;
        for ji=2:size(um,1)
            for jj=2:size(um,2)
                ua=0.5*((u2m(ji,jj)-um(ji,jj)*um(ji,jj))+(u2m(ji-1,jj)-um(ji-1,jj)*um(ji-1,jj)));
                va=0.5*((v2m(ji,jj)-vm(ji,jj)*vm(ji,jj))+(v2m(ji-1,jj)-vm(ji-1,jj)*vm(ji-1,jj)));
                
                eke_tot(ji,jj)=0.5*(ua+va);
                
            end
        end

%         % -- From ssh
%          [lon,lat,ssh,time_sec,time_txt,deptht]=rnemo_2d(bbio,bbio_dir,'gridT','ssh',year);
%          
          model.lon=lon;
          model.lat=lat;
%         
%         % ... Warning: Remove the equatorial band !!
%         eke_tot(1:size(ssh,1),1:size(ssh,2))=0.;
%         for it=1:size(ssh,3)
%             a=squeeze(ssh(:,:,it))-mean(ssh,3);   
%             a(find((model.lat > -2) & (model.lat < 2)))=NaN;
%             res_eke=eke(model.lon,model.lat,a.*100);
%             
%             eke_tot=res_eke+eke_tot;
%         end
        
        model.data=eke_tot;  
    case 'MLD'
        % MLD ... Delta sigma = 0.03 kg.m-3
        gridtype='gridT';
        
        
        [status,liste] = unix(['ls -al ',bbio_dir,'*',bbio,'*y',num2str(year),'*',gridtype,'.nc ']);
        % Indice de reperage pour les dates de chaque fichier
        txt=[bbio,'_'];
        inddate=findstr(liste,txt);
        
        
        model=struct;
        model.data(1:529,1:486,1:size(inddate,2))=NaN;
        
        for id=1:size(inddate,2)
            clear T_raw S_raw
            datet1=liste(inddate(id)+size(txt,2):inddate(id)+size(txt,2)+10);
            month=str2num(datet1(7:8));
            day=str2num(datet1(10:11));
            yy=str2num(datet1(2:5));
            kind='T';  
            levmin=1;
            levmax=46;      
            [model.lon,model.lat,T_raw,time_sec,time_txt,deptht]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax,yy,month,day);
            kind='S';                
            [model.lon,model.lat,S_raw,time_sec,time_txt,deptht]=rnemo_3d(bbio,bbio_dir,gridtype,kind,levmin,levmax,yy,month,day);        
            
            sigma=sw_dens0(S_raw,T_raw);
            
                       
            for ii=1:size(model.lon,1)
                for ij=1:size(model.lon,2)
                    
                    deltasig=0.;
                    ik=2;
                    while ((deltasig < 0.03) && (ik < length(deptht)-1))
                        ik=ik+1;
                        deltasig=sigma(ik,ii,ij)-sigma(2,ii,ij);
                    end
                    if (ik < length(deptht)-1)
                        model.data(ii,ij,id)=deptht(ik);
                    end
                end
            end
            % time_txt -> matlab time
            model.time(id)=datenum(str2num(char(time_txt(:,2:5))),str2num(char(time_txt(:,7:8))),str2num(char(time_txt(:,10:11))));
        end
    
end

